Visium HD Mouse Brain Analysis with Seurat and BANKSY

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Install and Load required libraries

# packages required for Visium HD
#install.packages("hdf5r")
#install.packages("arrow")

library(Seurat)

## Loading required package: SeuratObject

## Loading required package: sp

## 'SeuratObject' was built under R 4.4.1 but the current version is
## 4.4.2; it is recomended that you reinstall 'SeuratObject' as the ABI
## for R may have changed

## 
## Attaching package: 'SeuratObject'

## The following objects are masked from 'package:base':
## 
##     intersect, t

library(SeuratData)

## ── Installed datasets ──────────────────────────────── SeuratData v0.2.2.9002 ──

## ✔ ifnb     3.0.0                        ✔ stxBrain 0.1.2

## ────────────────────────────────────── Key ─────────────────────────────────────

## ✔ Dataset loaded successfully
## ❯ Dataset built with a newer version of Seurat than installed
## ❓ Unknown version of Seurat installed

library(ggplot2)
library(patchwork)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Load required libraries
library(DropletUtils)

## Loading required package: SingleCellExperiment

## Loading required package: SummarizedExperiment

## Loading required package: MatrixGenerics

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following object is masked from 'package:dplyr':
## 
##     count

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: BiocGenerics

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union

## The following object is masked from 'package:SeuratObject':
## 
##     intersect

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:dplyr':
## 
##     first, rename

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice

## The following object is masked from 'package:sp':
## 
##     %over%

## Loading required package: GenomeInfoDb

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

## 
## Attaching package: 'SummarizedExperiment'

## The following object is masked from 'package:Seurat':
## 
##     Assays

## The following object is masked from 'package:SeuratObject':
## 
##     Assays

library(Seurat)
library(hdf5r)

## 
## Attaching package: 'hdf5r'

## The following object is masked from 'package:SummarizedExperiment':
## 
##     values

## The following object is masked from 'package:GenomicRanges':
## 
##     values

## The following object is masked from 'package:S4Vectors':
## 
##     values

library(arrow)

## 
## Attaching package: 'arrow'

## The following object is masked from 'package:BiocGenerics':
## 
##     type

## The following object is masked from 'package:utils':
## 
##     timestamp

Data Setup Download the Visium HD 3’ Mouse Brain dataset from 10x Genomics: https://www.10xgenomics.com/datasets/visium-hd-three-prime-mouse-brain-fresh-frozen

Store the data in the folder VisiumHD_mouse_brain.

# Set working directory to where data is stored
#setwd("~/Desktop/VisiumHD_mouse_brain")

# Attempt to load data (will fail if required files are missing)
#localdir <- "/brahms/lis/visium_hd/mouse/new_mousebrain/"

#object <- Load10X_Spatial(data.dir = localdir, bin.size = c(8, 16))

This will produce: Error in FUN(X[[i]], …): File not found Cause: ‘filtered_feature_bc_matrix.h5’ is not found in the specified directory. Solution: Unzip the Visium_HD_3prime_Mouse_Brain_binned_outputs.tar.gz file.

## Load Processed Binned Output again
setwd("~/Desktop/VisiumHD_mouse_brain/binned_outputs/square_008um")

data.dir <- "~/Desktop/VisiumHD_mouse_brain/binned_outputs/square_008um"

object <- Load10X_Spatial(
  data.dir = data.dir,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "Spatial",
  filter.matrix = TRUE
)

## Data Inspection & Visualization

# Confirm default assay
Assays(object)

## An object of class "SimpleAssays"
## Slot "data":
## List of length 1

DefaultAssay(object) <- "Spatial"

# Visualize counts
vln.plot <- VlnPlot(object, features = "nCount_Spatial", pt.size = 0) + 
  theme(axis.text = element_text(size = 4)) + NoLegend()

## Warning: Default search for "data" layer in "Spatial" assay yielded no results;
## utilizing "counts" layer instead.

## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Rasterizing points since number of points exceeds 100,000.
## To disable this behavior set `raster=FALSE`

count.plot <- SpatialFeaturePlot(object, features = "nCount_Spatial") + 
  theme(legend.position = "right")

vln.plot | count.plot

Normalize and Explore Expression

object <- NormalizeData(object)

## Normalizing layer: counts

# Visualize gene expression
p1 <- SpatialFeaturePlot(object, features = "Hpca") + 
  ggtitle("Hpca expression (8um)")
p1

Unsupervised Clustering with Sketching

object <- FindVariableFeatures(object)

## Finding variable features for layer counts

object <- ScaleData(object)

## Centering and scaling data matrix

# Downsample for scalability
object <- SketchData(
  object = object,
  ncells = 50000,
  method = "LeverageScore",
  sketched.assay = "sketch",
  features = VariableFeatures(object)
)

## Calcuating Leverage Score

## Attempting to cast layer counts to dgCMatrix

## Attempting to cast layer data to dgCMatrix

# Switch to sketch assay
DefaultAssay(object) <- "sketch"
object <- FindVariableFeatures(object)

## Finding variable features for layer counts

object <- ScaleData(object)

## Centering and scaling data matrix

object <- RunPCA(object, assay = "sketch", reduction.name = "pca.sketch")

## PC_ 1 
## Positive:  Calm2, Nrgn, Snap25, Ppp3ca, Ywhah, Cck, Atp1b1, Rtn1, Olfm1, Chn1 
##     Ppp3r1, Ndrg4, Uchl1, Hpca, Syt1, Rasgrp1, Atp2b1, Enc1, Fam131a, Ncdn 
##     Tubb2a, Slc17a7, Gria2, Snca, Scn1b, Nsf, Ptk2b, Camk2n1, Nrn1, Pcsk2 
## Negative:  Ptgds, Plp1, Ttr, Mbp, Mobp, Enpp2, Apod, Ecrg4, Dbi, Igfbp2 
##     Trf, Clic6, Mal, Kl, Cldn11, Cnp, Folr1, Mag, Gsn, Mog 
##     Cryab, Fth1, Cd63, Calml4, Car2, Vim, Sparc, Pltp, Ppp1r14a, Cd9 
## PC_ 2 
## Positive:  Ttr, Ecrg4, Clic6, Kl, Enpp2, Folr1, Igfbp2, Calml4, 2900040C04Rik, Lbp 
##     Sostdc1, F5, Cab39l, Hemk1, Bsg, Fxyd1, Ace, Rbp1, Slc4a2, Kcnj13 
##     Arl6ip1, Mia, Trpm3, Col9a3, Wfdc2, Msx1, Tbc1d9, Slc2a12, Sulf1, Gm19935 
## Negative:  Plp1, Mbp, Mobp, Mal, Cldn11, Mog, Mag, Cnp, Bcas1, Gatm 
##     Pllp, Olig1, Ermn, Tspan2, Ndrg1, Opalin, Agt, Ppp1r14a, Sparc, Ugt8a 
##     Tmem88b, Septin4, Dbndd2, Nkx6-2, Gpr37, Slc6a11, Apod, Evi2a, Resp18, Ctss 
## PC_ 3 
## Positive:  Plp1, Mbp, Mobp, Fth1, Mal, Cldn11, Mag, Cnp, Mog, Trf 
##     Dbndd2, Car2, Tspan2, Septin4, Cryab, Pcp4, Prkcd, Ppp1r14a, Pllp, Opalin 
##     Ugt8a, Sgk1, Gpr37, Ermn, Apod, Ptgds, Tmem88b, Bcas1, Gatm, Enpp2 
## Negative:  Mgp, Acta2, Tagln, Myl9, Crip1, Nnat, 6330403K07Rik, Cldn5, Igfbp7, Vtn 
##     Vim, Tm4sf1, Resp18, Ly6h, Cavin1, Ly6a, Fn1, Itm2a, Ly6c1, Rgs5 
##     Flna, Mfge8, Ahi1, Dcn, Ltbp4, Pltp, Slc2a1, Timp3, Ndn, Zcchc12 
## PC_ 4 
## Positive:  Nrgn, Plp1, Hpca, Fam131a, Mal, Cldn11, Itpka, Baiap2, Camk2n1, Rprml 
##     Mog, Mag, Snca, Atp1a1, Cnp, Enc1, Icam5, Cracdl, Chn1, Fth1 
##     Ppp3ca, Ddn, Nptxr, Olfm1, Trf, Vxn, Klhl2, Gria2, Cnksr2, Tmsb4x 
## Negative:  Prkcd, Pcp4, Pdp1, Ramp3, Tnnt1, Rora, Adarb1, Ptpn4, Ntng1, Pcp4l1 
##     Cplx1, Slc17a6, Uchl1, Plcb4, Zic1, Shox2, Ccdc136, Synpo2, Plekhg1, Grm1 
##     Gabra4, Kcnc2, Ndrg4, Tcf7l2, Rab37, Nefh, Sparc, Spock3, Amotl1, Nsmf 
## PC_ 5 
## Positive:  Resp18, Ndn, Nap1l5, 6330403K07Rik, Ly6h, Zcchc12, Scg2, Ahi1, Gaa, Sst 
##     Calb2, Hap1, Pnck, Gprasp2, Tmem130, Cartpt, Tac1, Ttr, Gap43, Fxyd6 
##     Gad1, Nnat, Clic6, Ngb, Gad2, Vat1l, Npy, Ecrg4, Tubb2a, Rtn1 
## Negative:  Acta2, Mgp, Tagln, Myl9, Crip1, Igfbp7, Cldn5, Vtn, Ly6a, Ly6c1 
##     Apod, Cavin1, Tpm2, Fn1, Rgs5, Vim, Tm4sf1, Ptgds, Ltbp4, Flna 
##     Slc2a1, Epas1, Prkcd, Atp1a2, Dcn, Timp3, Itm2a, Myh11, Nbl1, Ccn2

object <- FindNeighbors(object, assay = "sketch", reduction = "pca.sketch", dims = 1:50)

## Computing nearest neighbor graph

## Computing SNN

object <- FindClusters(object, cluster.name = "seurat_cluster.sketched", resolution = 3)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 50000
## Number of edges: 2206986
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7948
## Number of communities: 59
## Elapsed time: 7 seconds

## 11 singletons identified. 48 final clusters.

object <- RunUMAP(object, reduction = "pca.sketch", reduction.name = "umap.sketch", return.model = TRUE, dims = 1:50)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## UMAP will return its model

## 00:05:51 UMAP embedding parameters a = 0.9922 b = 1.112

## 00:05:51 Read 50000 rows and found 50 numeric columns

## 00:05:51 Using Annoy for neighbor search, n_neighbors = 30

## 00:05:51 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 00:05:54 Writing NN index file to temp file /var/folders/tf/ns_lb6tj5fdbr5rfnl4w8zvw0000gs/T//RtmpwHvu2M/file165b073d84d49
## 00:05:54 Searching Annoy index using 1 thread, search_k = 3000
## 00:06:03 Annoy recall = 100%
## 00:06:03 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 00:06:03 34 smooth knn distance failures
## 00:06:04 Initializing from normalized Laplacian + noise (using RSpectra)
## 00:06:05 Commencing optimization for 200 epochs, with 2458170 positive edges
## 00:06:05 Using rng type: pcg
## 00:06:18 Optimization finished

## Project Back to Full Resolution

object <- ProjectData(
  object = object,
  assay = "Spatial",
  full.reduction = "full.pca.sketch",
  sketched.assay = "sketch",
  sketched.reduction = "pca.sketch",
  umap.model = "umap.sketch",
  dims = 1:50,
  refdata = list(seurat_cluster.projected = "seurat_cluster.sketched")
)

## full.pca.sketch is not in the object. Data from all cells will be projected to pca.sketch
## Projecting cell embeddings
## Finding sketch neighbors
## Finding sketch weight matrix
## Transfering refdata from sketch
## Projection to sketch umap
## Running UMAP projection
## 00:08:19 Read 376419 rows
## 00:08:19 Processing block 1 of 1
## 00:08:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 00:08:20 2493 smooth knn distance failures
## 00:08:20 Initializing by weighted average of neighbor coordinates using 1 thread
## 00:08:21 Commencing optimization for 67 epochs, with 11275170 positive edges
## 00:08:38 Finished

## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from full.umap.sketch to fullumapsketch_

## Visualize Clustering

# Sketched UMAP
DefaultAssay(object) <- "sketch"
Idents(object) <- "seurat_cluster.sketched"
p1 <- DimPlot(object, reduction = "umap.sketch") + ggtitle("Sketched Clustering") + theme(legend.position = "bottom")

# Full data UMAP
DefaultAssay(object) <- "Spatial"
Idents(object) <- "seurat_cluster.projected"
p2 <- DimPlot(object, reduction = "full.umap.sketch") + ggtitle("Projected Clustering") + theme(legend.position = "bottom")

## Rasterizing points since number of points exceeds 100,000.
## To disable this behavior set `raster=FALSE`

p1 | p2

## Spatial Visualization
SpatialDimPlot(object, label = TRUE, repel = TRUE, label.size = 4)

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

## Warning: ggrepel: 18 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

# Highlight specific clusters
cells <- CellsByIdentities(object, idents = c(0, 4, 32, 34, 35))
p <- SpatialDimPlot(object,
                    cells.highlight = cells[setdiff(names(cells), "NA")],
                    cols.highlight = c("#FFFF00", "grey50"),
                    facet.highlight = TRUE, combine = TRUE
) + NoLegend()
p

Marker Gene Analysis

object_subset <- subset(object, cells = Cells(object[["Spatial"]]), downsample = 1000)

## Warning: Not validating Centroids objects
## Not validating Centroids objects

## Warning: Not validating FOV objects
## Not validating FOV objects
## Not validating FOV objects
## Not validating FOV objects
## Not validating FOV objects
## Not validating FOV objects

## Warning: Not validating Seurat objects

object_subset <- BuildClusterTree(object_subset, assay = "Spatial", reduction = "full.pca.sketch", reorder = TRUE)

## Reordering identity classes and rebuilding tree

markers <- FindAllMarkers(object_subset, assay = "Spatial", only.pos = TRUE)

## Calculating cluster 24

## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session

## Calculating cluster 22

## Calculating cluster 46

## Calculating cluster 26

## Calculating cluster 5

## Calculating cluster 38

## Calculating cluster 9

## Calculating cluster 17

## Calculating cluster 15

## Calculating cluster 29

## Calculating cluster 44

## Calculating cluster 20

## Calculating cluster 27

## Calculating cluster 40

## Calculating cluster 35

## Calculating cluster 31

## Calculating cluster 39

## Calculating cluster 12

## Calculating cluster 4

## Calculating cluster 3

## Calculating cluster 45

## Calculating cluster 8

## Calculating cluster 13

## Calculating cluster 2

## Calculating cluster 0

## Calculating cluster 11

## Calculating cluster 34

## Calculating cluster 7

## Calculating cluster 25

## Calculating cluster 47

## Calculating cluster 33

## Calculating cluster 19

## Calculating cluster 10

## Calculating cluster 1

## Calculating cluster 41

## Calculating cluster 43

## Calculating cluster 18

## Calculating cluster 36

## Calculating cluster 21

## Calculating cluster 30

## Calculating cluster 28

## Calculating cluster 16

## Calculating cluster 6

## Calculating cluster 14

## Calculating cluster 42

## Calculating cluster 23

## Calculating cluster 32

## Calculating cluster 37

top5 <- markers %>%
  group_by(cluster) %>%
  filter(avg_log2FC > 1) %>%
  slice_head(n = 5)

object_subset <- ScaleData(object_subset, assay = "Spatial", features = top5$gene)

## Centering and scaling data matrix

## Warning: Different features in new layer data than already exists for
## scale.data

DoHeatmap(object_subset, assay = "Spatial", features = top5$gene, size = 2.5) + theme(axis.text = element_text(size = 5.5)) + NoLegend()

6. Identifying spatially-defined tissue domains While the previous analyses consider each bin independently, spatial data enables cells to be defined not just by their neighborhood, but also by their broader spatial context.

In Singhal et al., the authors introduce BANKSY, Building Aggregates with a Neighborhood Kernel and Spatial Yardstick (BANKSY). BANKSY performs multiple tasks, but we find it particularly valuable for identifying and segmenting tissue domains. When performing clustering, BANKSY augments a spot’s expression pattern with both the mean and the gradient of gene expression levels in a spot’s broader neighborhood.

# Install BANKSY and SeuratWrappers if needed
#if (!requireNamespace("Banksy", quietly = TRUE)) {
  remotes::install_github("prabhakarlab/Banksy@devel")

## Using GitHub PAT from the git credential store.

## Skipping install of 'Banksy' from a github remote, the SHA1 (098ad3f3) has not changed since last install.
##   Use `force = TRUE` to force installation

#}
#if (!requireNamespace("SeuratWrappers", quietly = TRUE)) {
  remotes::install_github("satijalab/seurat-wrappers")

## Using GitHub PAT from the git credential store.

## Skipping install of 'SeuratWrappers' from a github remote, the SHA1 (a1eb0d8b) has not changed since last install.
##   Use `force = TRUE` to force installation

#}


library(Banksy)
library(SeuratWrappers)

object <- RunBanksy(object,
                    lambda = 0.8, verbose = TRUE,
                    assay = "Spatial", slot = "data", features = "variable",
                    k_geom = 50)

## Fetching data from slot data from assay Spatial

## Subsetting by features

## Warning in asMethod(object): sparse->dense coercion: allocating vector of size
## 5.6 GiB

## Computing neighbors...

## Spatial mode is kNN_median

## Parameters: k_geom=50

## Done
## Done

## Creating Banksy matrix

## Scaling BANKSY matrix. Do not call ScaleData on assay BANKSY

## Setting default assay to BANKSY

## Warning: Layer counts isn't present in the assay object; returning NULL

## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the SeuratWrappers package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

DefaultAssay(object) <- "BANKSY"
object <- RunPCA(object, assay = "BANKSY", reduction.name = "pca.banksy", features = rownames(object), npcs = 30)

## PC_ 1 
## Positive:  Plp1.m0, Mobp.m0, Apod.m0, Ptgds.m0, Mbp.m0, Mal.m0, Trf.m0, Bcas1.m0, Cldn11.m0, Cnp.m0 
##     Mog.m0, Dbi.m0, Ppp1r14a.m0, Mag.m0, Ndrg1.m0, Prr18.m0, Qdpr.m0, Gatm.m0, Tspan2.m0, Gfap.m0 
##     Pdlim2.m0, Tppp3.m0, Ermn.m0, Gsn.m0, Plaat3.m0, Cryab.m0, Opalin.m0, Pllp.m0, Tmem88b.m0, Ugt8a.m0 
## Negative:  Calm2.m0, Ppp3r1.m0, Chn1.m0, Ppp3ca.m0, Ywhah.m0, Aldoa.m0, Atp1b1.m0, Snap25.m0, Rtn1.m0, Atp6v0c.m0 
##     Nptn.m0, Syt1.m0, Gria2.m0, Olfm1.m0, Dynll1.m0, Nrgn.m0, Clstn1.m0, Cck.m0, Cx3cl1.m0, Fam131a.m0 
##     Enc1.m0, Ptk2b.m0, Nsf.m0, Dnm1.m0, Atp2b1.m0, Mapk1.m0, Pcsk2.m0, Stmn3.m0, Ndrg4.m0, Syp.m0 
## PC_ 2 
## Positive:  Prkcd.m0, Adarb1.m0, Rora.m0, Ptpn4.m0, Pcp4.m0, Pdp1.m0, Tnnt1.m0, Ramp3.m0, Ntng1.m0, Synpo2.m0 
##     Rab37.m0, Plekhg1.m0, Plcb4.m0, Patj.m0, Shox2.m0, Pcp4l1.m0, Ccdc136.m0, Amotl1.m0, Cplx1.m0, Rgs16.m0 
##     Gabra4.m0, Zic1.m0, Nefh.m0, Nrip3.m0, Grm1.m0, Tcf7l2.m0, Elmo1.m0, Ptpn3.m0, Kcnc2.m0, Lhfp.m0 
## Negative:  Ddn.m0, Nrgn.m0, Baiap2.m0, Icam5.m0, Gda.m0, Itpka.m0, Cracdl.m0, Foxg1.m0, Snca.m0, Nptxr.m0 
##     Ppp1r1a.m0, Kcnj4.m0, Rprml.m0, Fam131a.m0, Egr3.m0, Ly6h.m0, Mpped1.m0, Gria2.m0, Ngef.m0, Stx1a.m0 
##     Ldha.m0, Cnksr2.m0, Ccn3.m0, Slc30a3.m0, Arpp21.m0, Rab40b.m0, Rasgrf2.m0, Lmo3.m0, Kcnv1.m0, Dgkb.m0 
## PC_ 3 
## Positive:  Lbp.m0, Calml4.m0, Mia.m0, Dynlrb2.m0, Gm19935.m0, Pltp.m0, Vim.m0, Fxyd1.m0, Enpp2.m0, Rsph1.m0 
##     Ttr.m0, Foxj1.m0, Igfbp2.m0, Rbp1.m0, Cd63.m0, Mlf1.m0, Csrp2.m0, Msx1.m0, Rarres2.m0, Cfap126.m0 
##     Prelp.m0, Pcolce.m0, Plxnb2.m0, Dbi.m0, Ezr.m0, Tm4sf1.m0, Pifo.m0, Ifi27.m0, Folr1.m0, Ecrg4.m0 
## Negative:  Prkcd.m0, Nsmf.m0, Rora.m0, Ptpn4.m0, Adarb1.m0, Ntng1.m0, Pdp1.m0, Tnnt1.m0, Ramp3.m0, Rab37.m0 
##     Plekhg1.m0, Rgs16.m0, Gabra4.m0, Shox2.m0, Pcp4.m0, Synpo2.m0, Kcnc2.m0, Grm1.m0, Cplx1.m0, Zdhhc22.m0 
##     Patj.m0, Ptpn3.m0, Elmo1.m0, Gabrd.m0, Ccdc136.m0, Zmat4.m0, Abhd12b.m0, Plcb4.m0, Dbndd1.m0, Nrip3.m0 
## PC_ 4 
## Positive:  Resp18.m0, Zcchc12.m0, Ndn.m0, Nap1l5.m0, Ahi1.m0, Gaa.m0, Scg2.m0, 6330403K07Rik.m0, Baiap3.m0, Gprasp2.m0 
##     Hap1.m0, Pnck.m0, Peg3.m0, Ngb.m0, Pgrmc1.m0, Ly6h.m0, Tmem130.m0, Impact.m0, Bex1.m0, Sparc.m0 
##     Nnat.m0, Pcbd1.m0, AW551984.m0, Calb2.m0, Gap43.m0, Gpx3.m0, Nrsn2.m0, Fxyd6.m0, Timp2.m0, Tmem91.m0 
## Negative:  Fth1.m0, Plp1.m0, Dbndd2.m0, Cldn11.m0, Mal.m0, Cryab.m0, Cnp.m0, Mobp.m0, Mag.m0, Mog.m0 
##     Trf.m0, Mbp.m0, Septin4.m0, Ppp1r14a.m0, Bcas1.m0, Car2.m0, Qdpr.m0, Sgk1.m0, Tspan2.m0, Ermn.m0 
##     Gatm.m0, Ndrg1.m0, Pllp.m0, Tmem88b.m0, Apod.m0, Ugt8a.m0, Opalin.m0, Prr18.m0, Gpr37.m0, Slc44a1.m0 
## PC_ 5 
## Positive:  Ndn.m0, Nap1l5.m0, Tppp3.m0, Resp18.m0, Zcchc12.m0, Qdpr.m0, Ahi1.m0, Plp1.m0, 6330403K07Rik.m0, Ly6h.m0 
##     Cnp.m0, Mal.m0, Mog.m0, Gaa.m0, Impact.m0, Gap43.m0, Cldn11.m0, Scg2.m0, Mag.m0, Tmem130.m0 
##     Dbndd2.m0, Gprasp2.m0, Mbp.m0, Bex1.m0, Mobp.m0, Stmn1.m0, Tspan2.m0, Pnck.m0, Trf.m0, Peg3.m0 
## Negative:  Igfbp2.m0, Nsmf.m0, Mgp.m0, Adcy1.m0, Ddn.m0, Fn1.m0, Cldn5.m0, Camk2n1.m0, Cavin1.m0, Slc6a20a.m0 
##     Rgs4.m0, Igfbp7.m0, Col1a2.m0, Vtn.m0, Colec12.m0, Rora.m0, Gabrd.m0, Myl9.m0, Acta2.m0, Tagln.m0 
##     Fmod.m0, Mef2c.m0, Pcolce.m0, Rbp1.m0, Cfh.m0, Slc22a6.m0, Slc6a13.m0, Krt12.m0, Tm4sf1.m0, Rgs5.m0

## Warning: Key 'PC_' taken, using 'pcabanksy_' instead

object <- FindNeighbors(object, reduction = "pca.banksy", dims = 1:30)

## Computing nearest neighbor graph

## Computing SNN

object <- FindClusters(object, cluster.name = "banksy_cluster", resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 376419
## Number of edges: 9391727
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9473
## Number of communities: 31
## Elapsed time: 50 seconds

## Visualize BANKSY Results

Idents(object) <- "banksy_cluster"
p <- SpatialDimPlot(object, group.by = "banksy_cluster", label = TRUE, repel = TRUE, label.size = 4)

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

banksy_cells <- CellsByIdentities(object)
p <- SpatialDimPlot(object,
                    cells.highlight = banksy_cells[setdiff(names(banksy_cells), "NA")],
                    cols.highlight = c("#FFFF00", "grey50"),
                    facet.highlight = TRUE, combine = TRUE
) + NoLegend()
p

Visium HD Mouse Brain Analysis with Seurat and BANKSY

Copy form Seurat tutorial

2025-08-10

R Markdown